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In this paper, we study Klein tunneling in random media. To this purpose, we simulate the 
propagation of a relativistic Gaussian wavepacket through a graphene sample with randomly dis- 
tributed potential barriers (impurities). The simulations, based on a relativistic quantum lattice 
Boltzmann method, permit to compute the transmission coefficient across the sample, thereby pro- 
viding an estimate for the conductivity as a function of impurity concentration and strength of 
the potentials. It is found that the conductivity loss due to impurities is significantly higher for 
wave-packets of massive particles, as compared to massless ones. A general expression for the loss of 
conductivity as a function of the impurity percentage is presented and successfully compared with 
the Kozeny- Carman law for disordered media in classical fluid dynamics. 
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INTRODUCTION 



As opposed to classical quantum mechanics where elec- 
trons tunneling into a barrier are exponentially damped, 
relativistic scattering was shown by Klein in 1929^ to fol- 
low a very unexpected behavior: If the potential is of the 
order of the electron mass or higher the barrier becomes 
virtually transparent to the electrons. This is called the 
Klein paradox. Experimental realizations were not avail- 
able until the recent discovery of graphene^^^. This ma- 
terial has revealed a series of amazing properties, such 
as ultra-high electrical conductivity, ultra-low shear vis- 
cosity to entropy ratio, as well as exceptional structural 
strength, as combined with mechanical flexibility and op- 
tical transparency. Many of these fascinating properties 
are due to the fact that, consisting of literally one single 
carbon monolayer, graphene represents the first instance 
ever of a truly two-dimensional material (the "ultimate 
fiatland"^). Moreover, due to the special symmetries of 
the honeycomb lattice, electrons in graphene are shown 
to behave like an effective Dirac fluid of massless chiral 
quasi-particles, propagating at a Fermi speed of about 
vf ^ c/300 ~ 10^ m/s. This configures graphene as an 
unique, slow-relativist ic electronic fluid, where many un- 
expected quantum-electrodynamic phenomena can take 
place, ^. For instance, since electrons are about 300 times 
slower than photons, their mutual interaction is propor- 
tionately enhanced, leading to an effective fine-structure 
constant agr = e'^/hvp ^ l- As a result of such strong 
interactions, it has been recently proposed that this pecu- 
liar 2D graphene electron gas should be characterized by 
an exceptionally low viscosity /entropy ratio (near-perfect 
fluid), coming close to the famous AdS-CFT lower bound 
conjectured for quantum-chromodynamic fluids, such as 
quark-gluon plasmas^. This spawns the exciting prospect 
of observing electronic pre-turbulence in graphene sam- 



ples, as first pointed out in Ref.-^ and confirmed by recent 
numerical simulations^. 

The zero-mass of electronic excitations in graphene 
may have other spectacular consequences. For instance, 
it has been recently pointed out^ that graphene could 
offer the first experimental of the so-called Klein para- 
dox, i.e. the capability of quantum wavefunctions to un- 
dergo zero reflection from a potential barrier much higher 
than the energy of the wavefunction itself. This prop- 
erty, which relies exclusively upon the spinorial nature of 
the Dirac wavefunction, stands in stark contrast with the 
corresponding non-relativistic behavior, which predicts 
an exponential decay of the transmission coefficient with 
the difference Vq — Vq being the height of the barrier 
and E the wavefunction energy. Based on an analytical 
solution of the scattering problem for a monochromatic 
plane wave, the authors were able to show that, depend- 
ing on a series of geometrical and energy parameters, 
special angles of incidence (resonant angles) provide lit- 
erally zero reflectivity: the plane wave goes completely 
across the barrier. Besides its intellectual charm, such 
property is of great practical interest for the study of 
electronic transport in graphene^' and it is expected 
to play an important role in the understanding of the 
minimum conductivity of graphene^^. Furthermore, the 
electronic spectrum of graphene can change depending on 
the substrate, for instance on SiC the energy spectrum 
presents a gap of width 2mvp^ which makes it possible to 
model the electric transport by using the massive Dirac 
equation^^^^^. Therefore, it is also interesting to study 
the Klein tunneling in a random media for this kind of 
gaped-samples (massive fermions case). 

On the other hand, due to the fact that, under suit- 
able conditions^, electronic excitations in graphene be- 
have as an effective relativistic Dirac fluid, in the pres- 
ence of a random media, transport laws similar to the 
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FIG. 1. The transmission coefficient of a Gaussian 
wavepacket, as computed with the analytical convolution, 
Eq. ()15|) . as a function of the incidence angle for a jD — 
0.15,0.31,0.46,0.92,1.85. The blue line corresponds to the 
unfiltered case, a ^ oo, corresponding to a plane wave. 

ones ruling fluid motion in diluted porous media, may be 
expected to apply. We refer here, e.g. to the Carman- 
Kozeny law^^^^^, which relates the permeability of a 
porous medium (conductivity of a graphene sample) to 
the solid concentration (impurity density). 

The paper is organized as follows: first, we introduce 
a brief description of the quantum lattice Boltzmann 
(QLB) method^^; second, we study the case of Klein tun- 
nehng of a Gaussian wave packet through a rectangular 
potential barrier. Subsequently, we present numerical so- 
lutions of the Dirac equation in the presence of random 
impurities, thereby providing an estimate for the effects 
of the impurity concentration on the conductivity of the 
graphene sample, for both cases, massless and massive 
Dirac fermions. The simulations are performed using a 
QLB model, which is also introduced as a new tool to 
study transport phenomena in graphene. Finally, we dis- 
cuss and summarize the results. 



II. THE QUANTUM LATTICE BOLTZMANN 
METHOD 

The quantum lattice Boltzmann (QLB) method^^ is a 
quantum- kinetic technique that was originally devised for 
non-relativistic quantum problems and recently shown 
to provide a second-order accurate solver for relativistic 
wave scattering and propagation^'^. Since the method 
is relatively new in the relativistic context, for the sake 
of self-containedness, we revisit here its main technical 
aspects. For full details, see our recent work^^^^^. 

The quantum lattice Boltzmann equation was initially 
derived from a formal parallel between the kinetic lat- 
tice Boltzmann equation and the single-particle Dirac 
equation— i^i^. For our purpose, it proves expedient to 
transform the standard form of the Dirac equation into 
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FIG. 2. The transmission coefficient of a Gaussian wavepacket 
as a function of the incidence angle for a = 24, 48 and 96 
(in lattice units), corresponding to cr/D = 0.46, 0.92, 1.85, as 
computed via convolution (solid line) and by QLB simulations 
(line with dots). 

the Majorana form, in which all matrices are real^, 

\dt + c(-a"S, + ^dy - a'd,) + iuJcay - igl] = 0, 

(1) 

This form is obtained by multiplying the standard Dirac 
equation on the left and right by the involution matrix 
U = 2~^/^(a^ + /3). In the above, c is the light speed, H is 
the reduced Planck's constant, / is the identity operator, 
and Uc = mc^/h is the Compton frequency for a particle 
of mass m. The wavefunction is a complex four-spinor, 
and a and /3 are the standard Dirac matrices. The last 
term couples the wavefunction to an applied scalar poten- 
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FIG. 3. Snapshots of the wavepacket density at various in- 
stants, t = 0, 420, 1050 (lattice units), for the case = (left) 
and = 27r/9 (middle), and (j) = 7r/3 (right) for cr / D = 1.85. 
In the middle, as one can see, after significant distortion in the 
intermediate stage of the evolution, the wavepacket manages 
to be transmitted across the barrier to a substantial extent 
(T = 0.76). On the other hand, at the right, the packet is 
mostly bounced-back by the barrier, with transmission coeffi- 
cient as low as T = 0.13. For visualization purposes, the color 
bar scale has been modified independently for each figure. 



tial V{x^ z) via the coefficient g = qV/h^ where q is the 
electric charge^^. Note that since the spin states mix-up 
during propagation (spinning particles), there is no basis 
in which all three matrices a^, a^, are simultaneously 
diagonal. 

Let us consider a one-dimensional version of the Dirac 
equation. In particular, let Z be a unitary matrix, diag- 
onalizing the streaming matrix —a^: 



1 

V2 



1 
1 

Vi 



1^ 
-1 

1 

1 0/ 



(2) 



Applying the matrix Z to Eq. (|T|), the streaming matrix 
along z is diagonalized and the collision matrix is also 
transformed accordingly 

[dt + cZ-\-a^)Zd, 

+ Z-\-ca''da: + cl3dy + iiOca^ - igI)Z]Z-^i^ = 0. 

(3) 

Neglecting any dependence of ip on the x and y co- 
ordinates, Eq. ([3]) may be written as a pair of one- 
dimensional Dirac equations 



dtdi^2 - cdzdi^2 = -0JcU2,i + igdi^2, 



(4) 



for the variables (1^1,(^2) and {u2,di) that represent the 
rotated wavefunction Z~'^ip = (ixi, 2x2, <^2)"^- The 
components u and d propagate up and down the z axis 
respectively, and the subscripts indicate the spin up (1) 
and spin down (2) states, respectively. The system of Eq. 
(jH) may be treated as a Boltzmann equation for a pair of 
complex distribution functions ui^2 and di^2^^. Equation 
^ may thus be discretized using the same approach as 
in lattice Boltzmann method, i.e. by integrating along 
the characteristic light-cones dz = :^cdt. 

The resulting system of algebraic equations reads as 
follows 



1 ^ 1 ^ 

Si, 2 - ^1,2 = 2^(^2,1 + C^2,l) + ^W{U1,2 + Si, 2), 

5i,2 - di^2 = -^fh{u2,i + ^2,1) + ^ig{di,2 + ^1,2), 



(5) 



where the hat superscript (^) indicates that the wave- 
function is evaluated at the end-point of the correspond- 
ing streaming step, namely 



Si, 2 = ^l,2(^ + t + At), Ui^2 = ^1,2(^, t) 

^1,2 = di^2{z - Az,t + At), di^2 = di^2{z,t). 



(6) 



The dimensionless Compton frequency is m = UcAt, and 
the dimensionless scalar potential is ^ = g{z,t)At. 

The pair of equations ([5]) can be solved algebraically, 
delivering explicit expressions for Si ^2 and ^1^2- 



Si, 2 = aui^2 + ^<^2,i, 
(ii,2 = adi^2 — bu2,i, 

where the coefficients a and b are 
1-1^/4 m 



(7) 



a 



ig' 



ig' 



m 



These coefficients satisfy |ap + |6p = 1, so that the right 
hand side of Eq. ([7]) corresponds to multiplying the ro- 
tated wavefunction = {ui,U2,di,d2)^ by the uni- 
tary collision matrix 



/a b\ 
a b 
-b a 
\-b ay 



(8) 



The streaming step propagates 1^1,2 upwards and di^2 
downwards, along the light cones given by Az = ±cAt. 
Note that this unitary operation is numerically exact, 
without round-off error, because the distribution func- 
tion is integrally transferred from the source to the desti- 
nation site, and no fractional transport is involved. Since 
both streaming and collisions step are unitary, the overall 
QLB scheme evolves the discrete wavefunction through a 
sequence of unitary operations for any value of the dis- 
crete time step At. In addition, since streaming pro- 
ceeds upwind only (no centered spatial differences) along 
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the discrete hght-cones associated with each component 
"^i , the QLB dispersion relation is automatically free 
from fermion-doubling,^^. This, together with the ex- 
cellent efficiency of the method, especially on parallel 
computers^^, should make QLB a potentially appealing 
candidate for computational studies of electron transport 
in graphene. 

The scheme extends to multiple dimensions through 
an operator splitting technique. Within this method, 
the three-dimensional Dirac equation splits into the sum 
of three one-dimensional equations, each involving spa- 
tial derivatives along one single direction. Each of the 
three stages representing evolution by a timestep dt is 
accomplished by rotating i/j to diagonalise the relevant 
streaming matrix, taking one timestep of the existing 
one-dimensional QLB scheme described above, and ro- 
tating back to its original basis. The algorithm is thus 
composed of the following three steps: 1) Rotate tjj with 
collide with X~^QX^ stream along x, rotate back 
with X; 2) Rotate with Y'^, collide with Y-^QY, 
stream along rotate back with Y; 3) Rotate i/j with 
collide with Z~^QZ, stream along rotate back 
with Z. This form emphasizes the symmetry between 
the three steps, but since the streaming matrix along y 
is already diagonal in the Major ana form, Y = I is the 
identity matrix. The matrix X reads as follows: 



X 



1 

71 



/-I 1 
10-1 
10 10 

V 1 1 



(9) 
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FIG. 4. Sketch of the domain setting used in our simula- 
tions of the propagation of a Gaussian wave packet through 
a porous medium. 



The rotated matrices X~^QX and Z~^QZ have the same 
sign pattern as Q, but Q does not. 

Summarizing, QLB provides a unitary, explicit algo- 
rithm for quantum wavefunctions in which information 
propagates along classical trajectories represented by a 
sequence of three one-dimensional light-cones, thereby 
avoiding any mixing of the spinorial components dur- 
ing the streaming step. Although detailed comparisons 
with other techniques remain to be developed, there are 
reasons to believe that such simplification may result in 
enhanced computational efficiency, especially with par- 
allel computers in mind. Finally, we wish to point out 
that the same algorithm describes both relativistic and 
non-relativist ic quantum wavepackets, depending on the 
value of the mass m and the characteristic strength of 
the potential energy. 



III. RELATIVISTIC GAUSSIAN 
WAVEPACKETS 



and the Z matrix is given in Eq. (|2]) above. 

The collision term splits into three parts, each of which 
is combined with the corresponding streaming step. The 
collision matrix thus coincides, up to a unitary transfor- 
mation, with the collision matrix for the one-dimensional 
QLB scheme, with a timestep ^dt (see Ref.^^). In par- 
ticular, Q is given by 



/a -b\ 
a b 
-b a 

\b a J 



(10) 



where the coefficients 
1 - ^^3/4 



a = 



b = 



ms 



are written in terms of the rescaled dimensionless Comp- 
ton and potential frequencies 



^3 = rrio 



ms = ^Ucdt, gs = ^gdt. 



The pattern of + and — signs in the b terms on the off- 
diagonal of Q follows the same pattern as the matrix. 



Since we are interested in apphcations on graphene, 
hereafter, our simulations will be performed in two spa- 
tial dimensions, (for more details see Ref.^^). The prop- 
agation of a plane wave through a rectangular potential 
barrier was discussed in Ref.^. However due to the fact 
that it only applies to monochromatic plane waves, i.e. 
infinitely extended states which may not necessarily be 
realized under all experimental conditions , it is of inter- 
est to explore to what extent are such results affected by 
the finite extent of the wavefunction. Here, for simplicity, 
we consider a Gaussian wavepacket of the form 



i^i{x,y) 



-f^^i{k^x+kyy)^ / = 1,2 (11) 



(4^a2)^/' 

where = + Ai = I /A, A2 = e*^/A with A 
^/A^~^~A^. The rectangular box potential of height Vq 
and width D is defined as follows: 



V{x) 



Vo, if < X < D, 
0, elsewhere. 



(12) 



Given the linearity of the Dirac equation and the fact 
that wavepackets are constituted by a Gaussian superpo- 
sition of plane waves, it is natural to express the trans- 
mission coefficient of a Gaussian wavepacket of size a 
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through the following convolution: 



k-k' 



CTk 



T{K,k'y)dKdk'y (13) 



where Sf = irk^^ with k\ = k^. ^ ky^ denotes the Fermi 
area, and G a Gaussian kernel of width cFk = 1/cr in 
mometum space. The function T(kx-, ky) is the transmis- 
sion coefficient of a plane wave with vector k = (kx^ky), 
which according to Ref. ^, can be calculated as T = 
1 - |rp with 

2ie''^{ss')-^ sin((7^D)(sin0 - ss' sin l9) 

^ ~ [e-'^-^ cos{(j) + 6>) + e'^-^ cos(0 - 0)] - 2i sm{qxD) ' 

(14) 

being <p the incidence angle, = {E — Vq)'^ / Ji^ ^ 
= tan~^{ky/qx) the refraction angle, s = sign(£^), 
s' = sign(£^ — Vo)^ and E the Fermi energy. 

Since the transmission coefficient for a plane wave only 
depends on the wave number ky , and due to the fact that 
the X component of the wave vector experiences a perfect 
transmission, as a first-order approximation^ we perform 
the convolution in just one dimension, ky, that is: 

T,{k)= / G(- — -]T{k')dk' (15) 

J-kF \ J 

where we have defined k = ky. By setting k' = k -\- 
and expanding T{k + q) around q = to second order, 
Eq. (p!5]) delivers 



T,(k)r^T{k)^^T"{k) 



0{al) 



(16) 



where T" is the second derivative of T with respect to 
k. The above expression means that resonant peaks 
{T'{kr) = 0, T"{kr) < 0) are smoothed out whenever 
the filter width a/c, is sufficiently high, or, more precisely. 



\T''(k)\ 
2T(/c) 



. This smoothing is the effect of non-resonant 
wave numbers. Given that a = 1/cr/c, one could readily 
estimate the minimal width a above which the secondary 
resonant peak would no longer be seen by the Gaussian 
wavepacket. However, the asymptotic expansion given 
by Eq. ([T6|) fails to represent the actual transmission co- 
efficient of the Gaussian wavepacket near the secondary 
resonant peak, the reason being that, around that peak, 
a second order expansion is grossly inaccurate because 
cr^T '^1 and higher orders will be even less accurate. 
As a result, the convolution integral, Eq. (p!5]) . needs to 
be computed. 



A. Computing the convolution 

To gain a quantitative sense of the dependence of the 
transmission coefficient of the Gaussian wavepacket with 
the spatial spread a, we have numerically computed the 
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FIG. 5. Transmission coefficient as a function of time for the 
impurity potential set at V = 50, 100, and 285 meV while 
varying the impurity percentage C {C =0.1%, 0.5%, 1% and 
5%) for m = 0. 



convolution integral of Eq. (p!5]) . for the following val- 
ues a/D = 0.15,0.31,0.46,0.92,1.85, where D = 100 
nm is the width of the potential barrier. The param- 
eters are the same as in Ref.^, namely E = 0.08 eV, 
Vb = 0.2 eV and D = 100 nm. The results are shown 
in FigdJ From this figure, it is seen that, for = 0, 
T{kr) = T{kFCOs{(j)r)) gocs from 1 to 0.7348, slightly over 
a 25 percent reduction. The same figure also shows that 
around the secondary resonance (at (j) = 27r/9), narrow 
wavepackets with a/D < 0.46 feature T ~ 0.5, with no 
sign of the secondary resonant peak. On the other hand, 
the secondary peak is seen to re-emerge for a/D > 0.92, 
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i.e. when cr is of the order of 100 nm, coraparable with 
the barrier width. With a / D = 1.85, the secondary peak 
is recovered, but only to about 80 percent. Note that, for 
(j) = 7r/2, the transmission coefficient is not zero, which is 
a consequence of the approximation made to obtain Eq. 
(p!5|) from Eq. (p!3|) . However, as shown in Sec. IIIIB| 
the numerical simulation of the transmission coefficient 
using QLB, shows generally a pretty satisfactory agree- 
ment with the approximation Eq. (fT5]) . 

In order to use the plane-wave approximation, one 
needs to ensure that the condition a > D \^ fulfilled, 
which sounds pretty plausible. However, this condition 
is strongly dependent on the angle of incidence. In partic- 
ular, it is far more stringent for oblique than for head-on 
{(j) = 0) incidence. Indeed, for ^ = 0, a/I) ~ 0.5 yields 
a substantial T = 0.9 for ^ = 0, while at (j) = 27r/9, 
we obtain a mere T ~ 0.4. At cr/D ~ 2, perfect trans- 
mission, T = 1, is practically recovered at ^ = 0, while 
for ^ = 27r/9, T - 0.8, i.e. about 80% percent of fuh 
transmission. 

We conclude that, for head-on incidence ~ 0), the 
transmission coefficient of Gaussian packets is still similar 
to the one of plane waves, as soon their extent becomes 
comparable to the barrier width. On the other hand, the 
secondary resonance, at oblique incidence, is highly af- 
fected by the finite-size of the wavepacket, and full recov- 
ery of perfect transmission seems to require wavepacket 
extents significantly larger than the barrier width. 



B. Numerical simulations 

The analytical expression of Eq. ([15]) has been com- 
pared against direct numerical simulation of the Dirac 
equation, using the quantum lattice Boltzmann (QLB) 
method. In order to back-up the previous findings, we 
have computed full numerical solutions of the Dirac equa- 
tion using a quantum lattice Boltzmann solver. The 
simulations are performed on 1024^,2048^ and 4096^ 
grids, depending on the size of the Gaussian packet. 
Lattice units are chosen such that D = D/Ax = 52, 
dt = dx/vp = 1.92 X 10~^^ seconds and energy is nor- 
malized in units of h/dt. The physical parameters are 
taken from Refi, that is E = 0.080 eV, Vq = 0.200 eV 
and D = 100 nm. The following sequence of wavepack- 
ets spreading, a = 24, 48, 96 has been simulated, with 
D = 52, all in lattice units. The results of the QLB sim- 
ulations appear substantially in line with the prediction 
of the convolution integral, i.e. they clearly show the dis- 
appearance of the secondary peak for a/D < 0.46, and 
its progressive reappearance above this threshold (see 
Figl2]). Note that, different from the solution of the con- 
volution integral, Eq. (p!5]) . the transmission coefficient 
measured by the simulation is zero for ^ = 7r/2, as should 
be expected, but the appearance of the second resonant 
peak is still retained. 

In Fig. O we show typical snapshots of the wavepack- 
ets for the cases ^ = 0, 27r/9, and 7r/3, for a/D = 1.85. 



The snapshots clearly show that, in the case ^ = 0, the 
wavepacket crosses the barrier totally unperturbed, with 
literally no distortion at any stage of the evolution. In 
the case of oblique resonant propagation, the packet still 
manages to cross the barrier to a large extent, (T = 0.76), 
with significant distortions in the intermediate stages of 
the evolution, leaving 24 percent of the packet behind. 
Finally, in the case of oblique non-resonant propagation, 
^ = 7r/3, the packet is mostly bounced-back by the bar- 
rier, with a transmission coefficient as low as T = 0.13. 



IV. KLEIN PARADOX IN RANDOM MEDIA 

One of the major technological challenges in cur- 
rent graphene research is to manufacture larger sam- 
ples, above 10 microns, for practical use in engineering 
devices^^. As the sample size is increased, however, it 
becomes more and more difficult to secure the purity of 
the sample, i.e. avoid crystalline inclusions (impurities) 
which alter the local structure of the graphene honey- 
comb lattice. Such impurities are indeed known to sig- 
nificantly affect the macroscopic properties of the sam- 
ple, primarily its electrical conductivity. To gain insight 
into this problem, it is therefore of interest to investi- 
gate the propagation of relativistic wavepackets within a 
disordered sample. 

The conductivity of two-dimensional massless fermions 
in disordered media has made the object of intense stud- 
ies in the literature, The contribution of the present 
work to this subject relates to the following three di- 
rections, i) Investigate the Klein-Paradox for the case 
of Gaussian wave-packets rather than plane waves, both 
for single barriers and disordered samples, ii) Discuss 
the viability of semi-classical descriptions of electrons 
excitations in disordered media, based on quantitative 
analogies with flows in porous media, iii) Expose the 
quantum lattice Boltzmann method as a new compu- 
tational tool for electron transport in graphene, which 
might bear a special interest for prospective implementa- 
tions on parallel computers. Notwithstanding points i-iii) 
above, we wish to point out that, being our solution based 
on the single-particle Dirac equation (no many-body ef- 
fects), any conclusion on transport phenomena in actual 
graphene samples must be taken with great caution. We 
also wish to remark that the Klein tunneling is expected 
to be relatively mild in the present set up, for two reasons. 
First, because the Gaussian wavepacket always includes 
non-resonant frequencies suffering partial reflection; sec- 
ond, because, being the wavepacket wider than the ob- 
stacle size (see below), it can split and turn around the 
obstacle like a classical fluid, hence be partially transmit- 
ted, without any quantum tunneling through the barrier. 

To analyze these transport phenomena, we simulate 
the propagation of a relativistic Gaussian wavepacket 
through a two-dimensional domain composed of three 
regions: an inlet region, where the wave packet is po- 
sitioned at the initial time t = 0; the impurity region. 



7 



1 






— C=0% 






— C=0.1% 


0.8 






— C=0.5% 








— C=1% 


0.6 






— C=5% 


0.4 






V = 50 meV 


0.2 


















1000 2000 3000 4000 5000 6000 
Time step 




C 



1 

0.8 
0.6 
0.4 
0.2 

0- 

1000 



— C=0% 

— C=0.1% 

— C=0.5% 

— C=1% 

— C=5% 



V = 100 meV 



2000 



3000 4000 
Time step 



5000 



6000 




— C=0% 

— C=0.1% 

— C=0.5% 

— C=1% 

— C=5% 



V = 285 meV 



1000 



2000 



3000 4000 
Time step 



5000 



6000 



FIG. 6. Momentum transmission coefficient Tjz as a function 
of time for tfie impurity potent iai set at V — 50, 100, and 
285 meV wfiiie varying tfie impurity percentage C {C =0.1%, 
0.5%, 1% and 5%) for m = 0. 



i.e. the central part of the domain where randomly dis- 
tributed barriers (impurities) are located; and the outlet 
region, which is the final region, where measurements of 
the transmitted wave packet are taken. Due to the large 
effective fine structure constant in graphene, we will ne- 
glect in our study the Coulomb interaction between carri- 
ers. The impurity concentration is given by C = N(P /A^ 
where N is the number of square obstacles of cross sec- 
tion (P, distributed over an area A = Ly x L^. For the 
present simulations d = S (larger than the typical lat- 
tice distante of graphene) and C is varied in the range 
0.001 ^ 0.05. In Fig. SI the computational domain is 



FIG. 7. Time at which 90% of the wave packet has been 
transmitted, to. 9, as a function of the impurity percentage 
for fixed values of V and m — 0. The potential barriers are 
as follows: V = 50, 100, 200 and 285 meV. The impurity 
percentage values are C =0.1%, 0.5%, 1% and 5%. 
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FIG. 8. Wave packet density p at times = 0, 900, 1500, and 
1800 (lattice units) for the simulation performed with impu- 
rity percentage C = 0.5% and V = 50 meV. 



sketched, periodic boundary conditions are imposed at 
top and bottom boundaries, while a bounce-back condi- 
tion is enforced at the inlet , and an open boundary con- 
dition is imposed at the outlet (so that the transmitted 
wave packet is not reflected back) . We use a square lattice 
of size 2048 x 512 cells, such that the regions [0, 512) x 512, 
[512, 1536) X 512, and [1536, 2048] x 512 correspond to the 
inlet, impurity, and outlet regions, respectively. The cell 
size is chosen to be Ax = 0.96nm, and the spreading 
of the initial Gaussian wave packet a = 48 (in lattice 
units), leading to a Fermi energy Ep = 0.117 (80meV in 
physical units). In our study, we use two values for the 
mass of the particles, m = (ungaped graphene) and 
m = 0.1 (gaped graphene), and vary the impurity po- 
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tential and the concentration. Five barrier heights are 
considered, namely V = 25,50,100,200,285 meV. Note 
that, while the first two lie below Ep^ hence can be over- 
come classically, the others can only be traversed head-on 
via quantum tunnelling. It should be further observed, 
though, that since the wavepacket is wider than the sin- 
gle impurity, i.e. a > even in the case Ep < V -> 
the wavepacket can split and turn around the obstacle 
like a "classical" fluid. Our results can be classified ac- 
cording to the energy of the particles, the potential of 
the barrier, and their mass as follows: weak potentials, 

V < Ep — mv\\ intermediate potentials, Ep — mv\ < 

V < Ef + mv'p] and strong potentials, V > E ^ mvp. 
The transmission coefficient T(t) is obtained by comput- 
ing T{t) = J^^^ ^ p{z^ t)dzdy^ where p is the wave 
packet density defined as p = \ui\'^ + |2i2p + + |<^2p, 
with i/j = (lii, 1^2, <ii, (^2)^ being the Dirac quadrispinor. 

A. Wave packet mass m — 

In this first set of simulations, we fix m = 0, and vary 
the impurity concentration, C, and the strength of the 
impurity potential, V. In Fig.O we fix the value of V and 
we compare T while varying the impurity percentage, in- 
cluding the reference value for the pure sample C = 0. 
From this figure, we observe that the wave packet takes 
longer to regroup for high impurity concentration and 
high impurity potential. This is a natural consequence 
of the randomness induced in the wavefunction by the 
disordered media. However, in all cases, the complete 
wave packet is reconstructed after some time, with no 
stagnant regions left behind. This can be related to the 
momentum loss due to the presence of the impurities, 
and therefore the motion of the wave packet experiences 
a corresponding slow-down. Note that, in order to re- 
cover the complete wave function, the simulations have 
been performed in a longer domain. Otherwise the right- 
moving wave-packet would leave the outlet region too 
early while the left-mover is still in the domain. In order 
to provide a measurement of momentum dissipation, i.e. 
the loss of conductivity due to impurities, we compute 
the momentum transmission coefficient as follows: 

Tj,{t)= I Jz{z,y,t) dzdy, (17) 

J Z>Zout 

where 

J,=i^^Aj+xf^Alxf, (18) 

is the z-component of the current density with Az the 
streaming matrix along z and = (i^i, 1^2, <^2)"^ the 
Dirac quadrispinor. 

In Fig. [6l we fix the value of V and compare Tjz^ while 
varying the impurity percentage. The subscript jz de- 
notes the transmission coefficient due to the z-component 
of the current density, J^. As a reference, we also plot 
Tjz{t) when the impurity percentage is set to C = 0. 
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FIG. 9. Transmission coefficient as a function of time for the 
impurity potential set at V — 50, 100, and 285 meV while 
varying the impurity percentage (C — 0.1%, 0.5%, 1% and 
5%) for m = 0.1. 

From Fig. [6] we can observe that, unlike the density, the 
momentum transmission coefficient does not saturate at 
unity (its value in the inlet region at the beginning of the 
simulation) , because momentum is irreversibly lost in the 
impurity region. Furthermore, as expected, the momen- 
tum loss increases with increasing impurity potential and 
concentration. 

As a characteristic quantity associated with the dy- 
namics of the transmission coefficient T, in Fig. [71 we 
report the escape time, to. 9, i.e. the time at which the 
transmission coefficient reaches 90%, (i.e. at 90% of the 
wave packet is transmitted through the obstacle region). 
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FIG. 10. Momentum transmission coefficient Tjz as a func- 
tion of time for tfie impurity potentiai set at 1/ = 50, 100, and 
285 meV wfiiie varying tfie impurity percentage {C = 0.1%, 
0.5%, 1% and 5%) for m = 0.1. 



As above, we plot to. 9 as a function of the impurity per- 
centage for two values of V. We notice that for high 
impurity concentration the Gaussian wave packet takes 
longer to cross the impurity barrier. The same effect oc- 
curs when the impurity potential is increased. At low 
impurity concentration, C = 0.001, the effect of the po- 
tential barrier is relatively minor, but, as the concentra- 
tion is increased, the escape time grows approximately 
linearly with the barrier voltage. 

In Fig. [HI we show some representative snapshots of the 
first 1800 time steps of the simulation, for impurity per- 
centage C = 0.5% and V = 50 meV. Here, we can see the 
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FIG. 11. Wave packet density p at times 0, 900, 1500, and 
1800 (lattice units) for the simulation performed with impu- 
rity percentage = 0.5% and V = 50 meV and with m = 0.1. 



way how the wave packet is scattered by the impurities, 
generating a plane front, as a result of the fragmentation 
of the wavefunction due to the random obstacles. 



B. Wave packet mass m = 0.1 

Next, we repeat the same simulations for the case 
of massive particles, with m = 0.1. Note that, since 
mVp/Ep = 0.83, the rest energy is a significant frac- 
tion of the kinetic energy, and therefore the wavefunction 
comes in the form of a superposition of two wavepackets, 
both moving at the Fermi speed, along opposite direc- 
tions, and mixing through the non-zero mass term. 

In Fig. [9l we fix the value of V and compare T, while 
varying the impurity concentration C. As a reference, we 
also plot T with C = 0. From the results, we observe that 
the wave packet takes longer to cross the impurity region 
than for the case of m = (the time it takes to reach a 
unit value of the transmission coefficient is longer). This 
is due the slow-down of the wavefunction as compared to 
the Fermi speed, because of the non-zero particle mass. 
Note the peak in the transmission coefficient, once the 
wave packet exits from the impurity region. This is due 
to the fact that Tj^ takes negative values in the late 
stage of the evolution, indicating the prevalence of the 
left-moving component of the wavepacket once the right- 
moving one has left the domain. 

We compute the momentum transmission coefficient 
using equations ([T7j) and (p!8|) . In Fig. [TOl we fix the 
value of V and compare Tjz while varying the impurity 
percentage. As a reference, we also plot Tjz{t) when 
the impurity percentage is set to zero. Note that, as 
expected, due to the inertia when the mass is increased. 
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the curve of the momentum transmission becomes wider 
than for the case of zero mass, reflecting the fact that 
the wave packet takes longer to move across the impurity 
region. In addition, the maximum momentum is smaller 
than for the case of zero mass, which indicates higher 
momentum losses. Thus, a non-zero mass of the (quasi)- 
particles, results in higher momentum losses. Also to be 
noted, are the negative values of Tj^ in the late stage 
of the evolution, indicating the presence of a left-moving 
component, most likely due to a spurious reflection at 
the outlet boundary. 

In Fig. [TTJ we show selected snapshots from the first 
1800 time steps of the simulation for impurity percentage 
C = 0.5% and F = 50 meV. From this figure, we observe 
that a portion of the wave packet gets "trapped", mov- 
ing at lower speed, within the impurity medium, while 
another portion manages to move out faster. 




C 



C. Momentum Transmission Coefficient Tjz 

In order to summarize the results obtained in the pre- 
vious sections, we inspect the maximum of the transmis- 
sion coefficient Tjz in Figs. [6l and [10) as a function of the 
impurity potential and concentration, for three different 
values of mass, m = 0, 0.05, 0.1 (see Fig. [12]). These data 
summarize the loss of momentum, hence resistivity, due 
to the random impurities, formally measured by 

??(C,y)=max(Tj,(C,F)) . (19) 

From these figures, we observe that at high impurity con- 
centration, C = 0.05, and a barrier V = 100 meV, the rel- 
ativistic wavepacket looses about 50% of its momentum, 
as compared the case of a pure sample {C = 0). At the 
same concentration, a massive wave packet with m = 0.1, 
would loose more than 80%, indicating a significant drop 
of transmissivity due to inertia. At low impurity level, 
C = 0.001, both massless and massive wave-packets show 
a mild reduction of transmittivity, below 10%. 
Let us now define the following "transmittance" : 

^{C,V) = -^ . (20) 

1 — 77 

This definition allows to draw a quantitative parallel with 
the concept of permeability of a classical fluid moving 
through a porous medium. That is, when the transmit- 
tance is unity, the conductivity goes formally to infinity, 
whereas zero transmittance connotes zero conductivity. 

Using Eq. ([2T]) . we have found that the numerical re- 
sults are satisfactorily fitted by the following analytical 
expression: 

nC,V) = A ^ +So , (21) 

where A, n, are fitting parameters, which depend on 
the strength of the potential and the mass of the par- 
ticles. In Fig. [T2I we report the results of the fitting 
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FIG. 12. Maximum value of Tjz as a function of the impurity 
percentage for each value of the impurity potential 1/ = 50 -r 
285. For three values of the mass, m = (top), 0.05 (middle), 
0.1 (bottom). 
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TABLE I. Set of parameters that has been obtained by fitting 
the numerical results for S using Ea. ()21|) . 

(solid line), showing good agreement with the numerical 
data. We have plotted 77 instead of S, in order to avoid 
the divergence at C = 0. The values of the parameters 
can be found in Table HI From this Table, we appreciate 
that the residual Sq, is zero when the mass is different 
from zero, which points to this minimum permeability 
(conductivity) as to a property of massless particles. On 
the other hand, massive particles show a closer adher- 
ence to the Kozeny-Carman law, in the context of clas- 
sical fluid dynamics^^^^^, where no residual conductivity 
is observed at C = 1. Also, note that for low potential 
barriers, the exponent is around n ~ 0.85, while for in- 
termediate and strong potentials it is near n ~ 1, i.e. 
the value it takes for classical fluid dynamics in a di- 
lute disordered medium. Thus, for strong potentials, the 
classical analogy shows satisfactory results, while for in- 
termediate and weak potentials, it presents deviations, 
typically of the order of 15%. Finally, we observe that 
the case m = shows a significantly higher transmission 
than the corresponding data with m > 0, which is due to 
the higher momentum losses in the impurity region. It 
appears plausible to interpret the non-negligible surplus 
of relativistic conductivity, especially for the three cases 
with Fermi energy Ep < as an indirect manifestation 
of Klein tunneling. 

V. CONCLUSIONS AND DISCUSSION 

In this paper we have performed a numerical study of 
a relativistic Gaussian wave-packet propagating through 
a disordered medium, which we modeled as a set of ran- 
domly located potential barriers. 

From the numerical results, we conclude that for high 
concentration of impurities, the wave packet presents 
higher losses in momentum. Furthermore, for a given im- 
purity concentration, by increasing the potential of each 
impurity, we also find a loss of momentum. Systems with 
massive excitations are also studied, which can be of rel- 
evance to the case of gaped graphene samples. A non- 



zero mass is found to produce higher losses of momen- 
tum in the impurity region. The actual numerical values 
show that at high impurity concentration, C = 0.05, the 
wavepacket looses more than half of its momentum with 
barriers of 100 meV and up to 85% with V = 285 meV. 
At low concentrations, C = 0.001, however, the losses are 
much milder, going from about 5 — 20%, for V = 100 to 
285 meV, respectively. 

These data can be regrouped into an analytical ex- 
pression, which bears a strong similarity with the per- 
meability of porous media, as a function of the porosity. 
We have estimated the value of the conductivity from 
the transmission coefficient and fitted it by using the 
Carman-Kozeny law for porous media, relating the per- 
meability with the concentration of impurities We have 
found that this analogy works pretty well for the massive 
case, which shows no residual conductivity and a scal- 
ing exponent pretty close to unity. On the other hand, 
the massless case shows a residual conductivity, which 
can possibly be related to the minimum conductivity of 
graphene. Moreover, for weak and intermediate poten- 
tial strengths, the exponent is not unity, corresponding 
to a fractional Kozeny-Carman law. On the other hand, 
for strong potentials, the exponent 1 is recovered to a 
good accuracy, bringing the results closer to the analogy 
with classical fluids^^. The applicability of this classi- 
cal analogy indicates that, at least for the parameter set 
investigated in this paper, quantum tunneling is not the 
dominant transport mechanism, as compared to the semi- 
classical dynamics of the wave-function, which can turn 
around the obstacles in a similar way as a classical fluid 
would do. The results of this paper are expected to be 
amenable to experimental validation. For this purpose, 
samples of graphene with local chemical doping could be 
used^'^^. In addition, for validating the results with mas- 
sive particles, a substrate of SiC will also be required, in 
order to generate the gap due to the presence of particle 
mass. Finally, as a byproduct, we have introduced a new 
tool to model electronic transport in graphene, namely 
the quantum lattice Boltzmann method (QLB). QLB 
shares a remarkable computational efficiency, especially 
on parallel computers, and easy handling of complex ge- 
ometries with its well-established classical LB counter- 
part. As a result, it is hoped and expected that the 
present model can make a contribution to the compu- 
tational study of transport phenomena in graphene and 
other physical systems governed by the Dirac equation. 
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